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Abstract 

Turbulence simulation codes can exploit the flute-like nature of plasma turbu- 
lence to reduce the effective number of degrees of freedom necessary to represent 
fluctuations. This is achieved by employing magnetic coordinates of which one 
is aligned along the magnetic field. In the most common implementation, posi- 
tions along the field lines are identified by the poloidal angle. 

In this work, an alternative approach is presented, in which the position 
along the field lines is identified by the toroidal angle. It will be shown that 
this approach has several advantages. Among these, periodicity in both angles 
is retained. 

This property allows moving to an equivalent representation in Fourier space 
with a reduced number of toroidal components. It will be shown how this duality 
can be exploited to transform conventional codes that use a spectral represen- 
tation on the magnetic surface into codes with a field-aligned coordinate. 

It is also shown that the new approach can be generalised to get rid of 
magnetic coordinates in the poloidal plane altogether, for a large class of models. 

Tests are carried out by comparing the new approach with the conven- 
tional approach employing a uniform grid, for a basic ion temperature gradient 
(ITG) turbulence model implemented by the two corresponding versions of the 
ETAI3D code. 

These tests uncover an unexpected property of the model, that localized large 
parallel gradients can intermittently appear in the turbulent regime. This leaves 
open the question whether this is a general property of plasma turbulence, which 
may lead one to reconsider some of the usual assumptions on micro-turbulence 
dynamics. 
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1. Introduction 



Microturbulence in tokamaks has the property that gradients parallel to 
the magnetic field are substantially smaller than gradients in the perpendicular 
direction. This is the consequence of low-frequency dynamics combined with 
a very efficient transport along field lines. On can assume for example that 
the typical evolution rate is of the order of the diamagnetic frequency lu* = 
(cT e / 'eB)kg I 'L, where L is a macroscopic (temperature or density) gradient scale 
length, and that parallel transport occurs at the rate given by free ion streaming 
(which can be considered typical of weakly collisional plasmas) 711 ~ Vthikn. 
Balancing these two rates one obtains 

k\\/k g ~ p s /L ~ p* (1) 

Codes have been developed that take advantage of this property to reduce the 
number of grid points to maximize code efficiency [2, [H, 11 3 ■ Moreover, using 
a coarse grid in the parallel direction reduces the time scale constraint given 
by the numerical stability condition, making the use of explicit schemes less 
penalising. The combination of these two factors gives to codes that are able 
to exploit this property an overall efficiency much superior than the one of 
conventional codes that employ a uniform grid and implicit schemes to deal 
with the parallel dynamics. 

This work explores ways to implement coarse gridding in the parallel direc- 
tion by making use of field- alignment. The initial motivation was the conversion 
of a certain class of conventional codes (spectral in the angles) to a more efficient 
representation. 

It turns out that there are just two practical ways to achieve this goal, which 
differ in the way one labels points along the field line. In the quasi-totality of 
codes that make use of field alignment, one employs a coordinate transforma- 
tion such that the parallel direction turns out to be effectively labelled by the 
poloidal angle. A second coordinate is needed to represent the fine structure 
of turbulent fields on a magnetic surface; this is effectively the toroidal angle. 
The alternative approach proposed in this work reverses the roles, by effectively 
using the toroidal coordinate to label the position of a point along the field lines, 
in a way that will be explained in detail in the following sections. The toroidal 
angle was previously employed as "parallel" coordinate in the "quasiballooning 
coordinates" construction of Dimits Q, which differs in important ways from 
the present work, as discussed later. 

It turns out that the approach pursued in this work has a number of advan- 
tages, a) It allows a simpler implementation of nonlinearities for a large class 
of models, b) It retains the original representation of functions in the poloidal 
plane, c) It can cope with X-point configurations including the open field line 
region, d) Functions are represented on a regular grid in real space. This allows 
moving to an equivalent representation in Fourier space with a reduced number 
of components. This duality can be exploited to transform conventional codes 
that use a spectral representation on the magnetic surface into codes with a 
field-aligned coordinate, e) Finally, as a further development, the new approach 
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Figure 1: Chart of electric potential fluctuations on a given magnetic surface 

can be generalised, for a large class of models, to get rid of magnetic coordinates 
in the poloidal plane altogether. Only the toroidal angle is needed to identify 
positions along the field lines to compute parallel gradients. 

This work is organised as follows. Section [2] reviews the main aspects of field 
alignment techniques. The new approach will then be derived in Sec. El giving 
first an intuitive description and then a mathematical justification. Section |4] 
will show how the spectral ITG code ETAI3D [H, Q can be converted to the 
new system with minimal programming. Section [5] is devoted to the comparison 
between the old and the new version of the code. These tests uncover an unex- 
pected property of the model. While gradients parallel to the magnetic fields are 
typically small, as expected as a consequence of the flute property fl| , localized 
sharp parallel gradients can appear in the turbulent regime. This is a property 
of the model (not of the method) which can have important consequences in 
terms of convergence. Generalisation to a system not using magnetic surfaces 
is described in the Appendix. 

2. Review 

The problem posed by the two-scale nature of plasma micro-turbulence is 
illustrated in Fig. [TJ 

This figure shows a chart of the turbulent electrostatic potential $ at a fixed 
radial position, that is, on a given magnetic surface, as a function of the toroidal 
angle ip (abscissa) and poloidal angle 8 (ordinate). Details of how this figure was 
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Figure 2: Temperature fluctuations spectrum at a given radial position 

obtained will be given later. One can clearly observe a direction along which 
structures tend to align along the lines 9 = ip/q(r) + const. , which correspond to 
the local magnetic field direction, q{r) being the safety factor at the given radial 
location r. It is clear that turbulent structures have generically small parallel 
gradients. Carrying out a Fourier transform to mode number space (n, m) one 
obtains the spectrum \$ m ,n\ 2 given in Fig. [2] 

As expected, most of the turbulent energy is concentrated along the lines 
m = nq. It is apparent how much computer resources are wasted by a repre- 
sentation of the turbulent fields that uses a grid with uniformly small spacing 
in the angles, as well as by its corresponding representation in mode number 
space. Clearly, a coordinate system that would exploit the property of small 
parallel gradients would constitute a substantial advantage. 

For the sake of simplicity, only magnetic configurations with circular concen- 
tric magnetic surfaces will be considered in this work. Usual radial and angular 
coordinates (r, 9, tp) are employed. Information on the local magnetic field di- 
rection is contained in the safety factor q(r). The derivative along the magnetic 
field is embodied in the parallel derivative operator: 



d 



V„ = — 



1 d 



dp q(r) 09 



(2) 



This operator is proportional to the actual parallel derivative operator by an 
r-dependent factor. The first implementation of field aligned coordinates in a 
plasma turbulence code can be traced back to thinking in the early Nineties Q 
that led to Ref. [lj ■ Consider the transformation: 



£ = tp - q(r)9 
s = 9 
p = r 



(3) 
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Derivatives with respect to the original variables are given by 

d d ,. Q d 

d-r = dp- q ' {r)6 dZ 

8_ _ 8_ 
dp d£ 

8 8 8 
89 = ~d~s~ q ^'d£ 

While the parallel derivative is given by 

1 8 



(4) 



q(r) 8s 



(5) 



Since parallel derivatives are expected to be small, a small number of points 
along s is needed to represent a function adequately. In the light also of Eq. [3J 
one can equivalently say that a small number of poloidal points is needed to 
represent a given function. One also observes that the poloidal derivative is 
given by a linear combination of a (slow) derivative along s and of a (fast) 
derivative along £. Since the latter is a derivative taken along ip, one can say 
that the transformation §5§ decomposes the poloidal derivative into a parallel 
component and a toroidal component. Since the parallel variable is coarse, in 
this representation all the information on the fine structure of turbulence is 
necessarily carried by the toroidal angle. Also, in actual code implementation, 
the derivative along s is usually dropped, approximating the poloidal derivative 
by 8/89 w q{r)8/8^. In linear theory, for a given toroidal mode n, one can 
Fourier-transform along ip (£), which gives 8/89 rs inq(r). One can see that 
transformation ((3]) is equivalent to the ballooning transformation Q, as already 
remarked in Q. 

This coordinate transformation, as it is, has some drawbacks. First, one 
notices that the new coordinate s is not periodic. Care must be taken, when 
setting boundary conditions at the end points of an s line, that the original 
double-periodicity of the given magnetic surface is enforced |3| . Non-compliance 
with this constraint would lead to spurious solutions and dubious results. The 
second problem is the consequence of the term, proportional to 9, appearing 
in the expression of the radial derivative Eq. 2) This term, familiar from the 
ballooning transformation, leads to mixed-derivatives of increasing weight as 
one moves away from = 0, when computing certain operators as the Lapla- 
cian. These terms are the consequence of ^-dependent non-diagonal metric 
coefficients in the new coordinate system. Although mathematically correct, 
they pose a numerical challenge since their numerical treatment can introduce 
artificial inhomogeneities in the poloidal direction even for system possessing 
poloidal symmetry. This problem is cured by the shifted- metric technique jj| 
as described below. A third problem with the old coordinate transformation 
is that it cannot deal with the separatrix, since, there, the safety factor be- 
comes infinite. This problem is automatically solved with the new coordinate 
transformation proposed in this work. 
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The shifted-metric approach introduced in [4{ consists in sectioning the 
toroidal manifold given by the magnetic surface into a number of parts, Ng, 
each having its own coordinate system differing one from the other by a shift 
in the origin. Here a slightly different version is given, which makes use of 
overlapping poloidal sectors. Consider the set of transformations 

t = <p-q{r){0-9 k ) 

Sk = e- e k (6) 

p = r 

where 6% = A9k with A8 — 2-k/Nq and k = 0, Ng — 1. A given sector is defined 
by the rectangle such that < <p < 2tt and 9 k — A9 < 9 < 9 k + A9. Derivatives 
with respect to the original variables are now given by 



_8_ _ d_ 

dip d£ 
d d d 
89 = !h~ q ^'dti 



(7) 



One observes that, for each poloidal sector k, the additional term that leads 
to metric distortion vanishes at 9 = 9 k - If now one uses exactly the same 
number of sectors as the number of points along 9 needed to achieve the required 
resolution, one realises that it is necessary to compute perpendicular operations 
only at points where no metric distortion occurs. The computation of parallel 
derivatives requires the knowledge of the functions to be derived at the end 
points sj; = ±A8 of an s k line. Whereas one could use collocation points at 
these ends, this would still leave the question of relating values of the function 
at these points with the values at the s k ±i = lines of the neighbouring k ± 1 
sectors. In practice, the solution adopted in which will be maintained here, 
is to assign values of a function only on the s k = lines, and, for each k, 
to compute the values at the end points s k = ±A9 by interpolation of the 
values at the s k ±i = lines of the neighbouring k ± 1 sectors. Notice that the 
accuracy needed to compute the values of a function at the interpolation points 
is automatically assured by the high resolution needed to describe a function on 
a given s k = line. 



3. New approach to field-aligned coordinates 

One starts with the observation that the net result of the shifted metric 
procedure is that functions are now defined on a regular grid in the former 
(if, 9) space, with the advantage that the discretisation is now coarser in 9 than 
it would be in a conventional meshing with uniformly small spacing in both 
directions. The coordinate 9 labels the position along the field lines. Also, this 
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Figure 3: Mesh on a magnetic surface with reduced resolution in 9 (old approach) 



mesh still has the good property of being periodic in both directions, which are 
effectively the original angles. The situation is described in Fig. [3] 

One also remarks that, as a consequence of the flute property, one can reduce 
the resolution in any chosen direction, provided that the information on the 
fine structure of the turbulent fields can be carried by the variation in any 
other direction. One then realises that there is just another alternative to this 
meshing, which preserves the good property of double periodicity. It is given by 
switching the roles of coarse/fine mesh between the poloidal/toroidal angles, as 
shown in Fig. 01 

This latter choice has a coarse toroidal grid. This implies that it is now 
the toroidal angle that must be used to describe variations along the field lines. 
This can be justified mathematically by sectioning the toroidal manifold in N v 
overlapping toroidal sectors, each with its own set of field-aligned coordinates 
given by the family of transformations 



q(r) 



Sk 

P = 



(8) 



where (pk = Atpk with Aip = 2ir/N ip and k — 0, N v — 1. Derivatives with respect 
to the original variables are now given by 
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Figure 4: Mesh on a magnetic surface with reduced resolution in tp (new approach) 



d _ d Id 
dip ds q(r) d£ ' 
d_ _ 8_ 

whereas the parallel derivative is now simply given by 

d 



(9) 



(10) 



On any given ipk = line, the derivatives are given by the same expressions as 
the original derivatives. Parallel derivatives are computed as before by interpo- 
lation at the end points of the ipk lines. 

This new system has additional advantages. The fact that the transforma- 
tion involves l/q(r) rather than q(r) allows one to take the limit of q — > oo 
without problems; this permits to treat a system whose magnetic equilibrium 
possesses a separatrix. Also, the poloidal derivative is simply given by a deriva- 
tive along £, while the radial derivative is a simple derivative along p (at (p = tpk). 
This permits to implement key operators like the Laplacian or perpendicular 
Poisson brackets in a manner strictly equivalent to the one of the initial coor- 
dinate system. This fact will be exploited in Sec. [5] to carry out a comparison 
between mathematically equivalent codes that employ the conventional and the 
field-aligned coordinate systems. 

As mentioned in the Introduction, the approach presented here has in com- 
mon with Ref. the use of the toroidal angle to keep track of the position 
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along the field lines. However, there are several important differences. Ref. [2J 
uses, for each radial point, an (optimal) rational approximation to \/q{r) in 
what would be Eq. [5] As a result, parallel derivatives expressed in terms of 
the new variables still contain a contribution from the derivative along the fast 
variable, albeit with a small coefficient. Moreover, points are collocated along 
the new almost-parallel direction, without use of the shifted-metric idea, and as 
a consequence one looses periodicity in one direction. Perpendicular operators 
like the Laplacian need to be computed with interpolation, unlike in the present 
work, where interpolation is necessary only to compute parallel derivatives. 

Finally, the approach presented here can be extended to avoid using magnetic 
coordinates in the poloidal plane, as sketched in the Appendix. 

4. Implementation of field-alignment in the ETAI3D code 

In order to test the new approach, one has first to implement it in some 
code that solves a three-dimensional micro-turbulence model with a sheared 
magnetic equilibrium as minimal features. It was found convenient to start from 
an existing code, ETAI3D, which solves a model of ITG turbulence described 
by the following equations: 

^w-2euj d ^ + Tw + Ti) + AV||«|| = V(D w Vw) , (11) 

^-vn -4r £ w d vn + AVii ($ + tw + TA =V(D Vun), (12) 

at " 11 " 11 

^Ti - 2 r e u> d ((r - 1)($ + tw) + (2T - 1)T,) + 

+A t (r - 1)V||«|| = -A lL r 1 / 2 |V|, | Ti + V(D T VTi) . (13) 

Here w is the ion guiding centre density, i>|| the ion parallel velocity and Tj 
the ion temperature. $ is the electrostatic potential, which is related to w by 
w = $— < <I> > — pJ.Vj^, where < . > denotes the flux surface average. The 
equations are implemented on a cylinder (r, 9, <p) where tp is the direction of the 
principal (would-be toroidal) magnetic field. 

The operators in Eqs. f| 1 ltil 3[) are: a) the total time derivative in the ExB 
velocity field ve, 4ff = {dt + ve • V)/ = d t j + [$,f] , where the Poisson 
bracket, in cylindrical coordinates, is [$,/] = d r ^^dgf — d r f~d$$>; b) the 
curvature operator cj^ = sin 8 d r + - cos 9 dg ; c) the parallel derivative operator 
V|| = d v + ^dg, with q the safety factor. 

The parameters of the model are the ion sound Larmor radius normalized 
to the minor radius = p s /a, the inverse aspect ratio e = a/R, A — e/p*, a 
constant T = 5/3, the temperature ratio r = Tio/T e Q, the small perpendicular 
transport coefficients D w , D v and Dt, set to damp the small scales, and 7l = 1 
is the coefficient of the parallel heat flux in Eq. (fT3"j) . modelled by a Landau- 
fluid closure of the Hammett-Perkins type [l(|, as appropriate in the weakly 
collisional, long mean- free-path limit. 
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Lengths are normalised to the minor radius a, times to the Bohm unit 
tBohm = a 2 / XBo hm, wit h XBohm = p s c s , the parallel velocity to the ion sound 
velocity c s = yjT e /m.i, the ion temperature to T e = T e o (assumed spatially 
constant) and the potential to T e /e. 

Note that this model, involving only three scalar fields, is substantially sim- 
pler than ITG models implemented in more modern codes that aim to reproduce 
a good agreement with observations. However it is adequate for testing purposes 
in the context of the present work. 

The fundamental variables evolved by the code are the Fourier components 
of the fields in the angular variables, truncated to given maximum poloidal 
and toroidal mode numbers, and defined at given points on a discrete radial 
grid. Derivatives are computed by second order finite differences in the radial 
direction and by direct multiplication by the corresponding wave numbers in 
the angular directions. Quadratic nonlinearities are computed with the pseudo- 
spectral method, by transforming first to real space, computing the products 
there, and coming back to Fourier space in the end. A de-aliasing technique as 



given in Ref. [11| is employed to ensure that the pseudo-spectral procedure to 
calculate nonlinearities is strictly equivalent to a direct convolution truncated at 
the maximum wave numbers. Time advancing makes use of a splitting technique 
where the reactive part of the model (nonlinearities and curvature) is advanced 
by a leap-frog scheme, and the diffusive and parallel operations by implicit 
schemes that keep the algorithm to overall second order accuracy. 

In this standard spectral representation, the parallel derivative operator ((2J 
is represented by n + m/q{r). As already shown in Fig. [2j this type of spatial 
discretisation is rather inefficient. 

Consider now the question of the implementation of the parallel derivative of 
a given function F(r, 9, y?) as obtained by the field-alignment transformation ®. 
At any given point in real space, one can use a second order centred finite 
difference (FD) expression of (fTOj) . computed with a suitable spacing As along 
the coordinate s, while holding £ and p constant. In the original coordinate 
system, the FD expression of the parallel derivative then reads: 



= Fjs + As) - Fjs - As) = 
11 2As 
F(r, 9 + Ay?/ g (r), y? + Ay?) - F(r, 9 - Ay?/q(r), y? - Ay?) 

2Ay> 



(14) 



where the fact that As = Ay> has been used. We note again that values of the 
function at the end points (y? ± Ay?, 9 ± Aip/q) along the field lines need to be 
computed by interpolation. For this, a good knowledge of the function along 
the poloidal direction (lines of constant y?) is necessary. That is, one must have 
in general A9 << Ay? to resolve the turbulence scale length in the poloidal 
direction, whereas Ay? needs not be very small. 

Since one works with a doubly periodic domain in the angles, interpolation 
is conveniently carried out by transforming to Fourier (mode number) space. 
The two representations have the same information content. This implies that 
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if a small number of toroidal points gives an adequate representation of a func- 
tion, the same (small) number of toroidal modes is needed, provided that the 
resolution in the poloidal direction is adequate. 

Consider the Fourier representation of a function F, written in terms of its 
components F mn . 

F(r, e,<p) = J2 F mn (r)e^ m9+n ^ (15) 

Once F is known on the grid nodes in real space, this expression can be used 
to obtain approximate values of the function at all points in space. The FD 
expression of the parallel derivative (fT4)l becomes 

1 — ' A(j9 

One then observes that the expression 

sin[(n + m/q(r))Aip] 
Aip 

is the representation of the parallel derivative operator in mode number space, 
in its approximate FD form. One notes that for small enough values of Aip, 
such that max[(n + m/q(r))Aip] « 1 this expression reduces to the usual one 
n + m/q(r). However, from the previous discussion, one expects that expression 
(fT7|) can be used also with substantially larger values of Aip, 

In order to see why it works, one must keep in mind that turbulent energy 
accumulates around regions where k\\ « 0. With the conventional spectral 
representation this occurs where n + m/q(r) » 0. One often refers to this 
occurrence as the " resonance condition" . With the representation given by (|17[) 
the condition is 

sin[(n + m/q(r))Aip] (18) 

which is satisfied by all combinations (n + m/q(r))Aip — kir with k integer, 
up to some maximum value. It is easy to see that when AO « Aip, the 
maximum allowed value of to, tt/ AO, is sufficiently large that the resonance 
condition occurs for several values of A:. As a consequence, for a given domain 
in (n, to) space, several bands satisfying ([T5]) are populated. The number of 
bands increases with the ratio Aip/ AO, while the total number of sites in (n, to) 
space that satisfy (jT5)) stays roughly constant. 

The situation is illustrated in Fig. [SJ which shows the spectrum of the same 
field as in Fig. [2] after a transformation to real space, reduction of the toroidal 
resolution by a factor 8, and transformed back to mode number space. This 
result can be seen as the effect of aliasing, which is beneficial in this instance, 
since aliased modes occupy formerly empty space. In this case one can say that 
aliasing acts as an efficient mean of data compression. 

When solving a plasma micro-turbulence model, the condition of small par- 
allel wave-number is produced dynamically. If the initial conditions are such 
that the modal composition satisfies the appropriate resonance condition, this 
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Figure 5: Temperature fluctuations spectrum in field-aligned coordinates 

property is preserved by the subsequent evolution. In particular, initial con- 
ditions with a spectrum like the one shown in Fig. [5] will continue to have a 
spectrum of this type if the parallel derivative operator is of the form given in 

(HZD- 

It has then been a simple matter to convert a code like ETAI3D to field 
aligned coordinates by replacing all the occurrences of the combination of n + 
m/q(r) by expression (|17[) . An important point to note is that the spectral 
nature of the code allows one to keep the Landau closure on the parallel heat 
flux in the original form proportional to V|| as given in Eq. I13L It is so also for 
the modified version with a coarser toroidal resolution, since | V|| | is replaced by 
the absolute value of (fT7]) . In this respect the new code is unique, since parallel 
derivatives in field-aligned codes are usually treated with finite differences in 
real space. This makes it difficult to deal with nonlocal dissipation operators 
proportional to |Vy |, which are then replaced by higher derivatives. It is impor- 
tant to keep in mind this peculiarity when examining the results of the tests of 
the new version of the code, presented in the following section. 

5. Tests 

Tests of the new version of ETAI3D were carried out at values of p* of 1/50 
and 1/100 in cylindrical mode (with the curvature operator u>d set to zero) and 
in toroidal mode, with various initial conditions. The tests consist in comparing 
the results by varying the number of toroidal points N v . All the tests were 
carried out with a simulation domain comprised between r = 0.5 and r = 1, 
and with a safety factor profile q = 2r. 

In general, similar results were obtained for all values of N v as low as 32. 
On observes that the behaviour at the lower value of N v = 16 can be very dif- 
ferent, at least transiently. The situation is illustrated in Fig. [6] The plots show 
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the amplitude of turbulent fluctuations of a cylindrical case with p* = 1 /50, ob- 
tained at Ng — 128 (fixed) and with N v of 128, 32, and 16. Initial conditions are 
those of an unstable temperature profile with the addition of small fluctuations 
that satisfy the small parallel gradient condition. The system initially evolves 
linearly, with the fluctuations growing exponentially, and eventually saturates 
in a turbulent state. 
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Figure 6: Comparison of the turbulent fluctuation amplitude of potential (solid) and temper- 
ature (dash), obtained at N v = 128 (blue lines), N v = 32 (green lines) and N v = 16 (red 
lines) 

The blue lines represent simulations at N v = 128. When compared with 
the former version of the code that does not employ field-aligned coordinates, 
the time traces are virtually undistinguishable (these are not shown here) . This 
is taken as the reference case. As the number of toroidal points is reduced to 
N v = 32 (green lines) one can note some difference in the nonlinear phase. 
This can be attributed to the reduced accuracy inherent in using a larger Atp 
to compute the parallel derivative. Upon a further reduction to N v — 16 (red 
lines), one observes that the traces strongly overshoot the reference case when 
entering the nonlinear phase, although eventually one observes signs of recovery. 

In order to understand why the code breaks down, consider that the following 
criterion must be satisfied at every point in the simulation domain 

k\\(r,6,ip)A<p < 1 (19) 

where fcy (r, 9, if) is a local measure of the parallel gradient. The following defi- 
nition is employed for any function F(r, 9, ip) 

fc,l(r,0,¥>)=V||J7||F|| (20) 



13 




k— parallel 

Figure 7: Histogram of fcii for a N v = 32 case 



where \ \F\ \ is the mean squared norm at radial location r, 



\F\ 
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Fig. [7] and Fig. |8] show the distribution function of fe|| of the iV^, = 32 and 
N v = 16 cases respectively, at the end of the simulation. One notices that 
the distribution of the N v = 16 case is broader, with wings at fairly high fen. 
Fig. [9] compares the maximum and the r.m.s. value of fen in the two cases. It is 
apparent that criterion (fl~9)) is satisfied in the N v — 32 case, since the highest 
value of fc|| is about 3 so that the maximum value of expression (|19p achieved at 
any point in the simulation domain is about 0.5. It is also clear that when N v 
is reduced to 16, the criterion is violated. When this happens, the discretised 
system is not able to adequately represent the continuum system. 

Several tests show that the linear phase is the easiest to simulate (the one 
that tolerates smaller N v ), while the most critical one is the transient phase 
after the linear one when the system adjusts to the saturated turbulence. In 
this phase, locally higher-than-expected values of fc|| can occur, although the 
typical value stays close to the inverse of the characteristic system length. 

It must be understood that the generation of locally strong parallel gradients 
depends on the model, not on the method. Indeed, analysis of older simulations 
carried out with the conventional uniform discretisation in the angles, also show 
this behaviour. An example drawn from data from the European code bench- 
marking on Cyclone test cases is given in Fig. [TDJ In this instance, 
values of k\\ as high as 40 (although with a minuscule probability) can occur. 
This would require a minimum of N v fs 240 to stay on the safe side. 

For any given plasma turbulence model it is obviously important to assess 
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Figure 8: Histogram of fcy for a Nip = 16 case 
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Figure 9: Maximum values (dashed lines) and r.m.s. values (solid lines) of ka as a function 
of time, for the two cases at Nip = 16 (green upper lines) and Nip = 32 (blue lower lines). 
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Figure 10: Histogram of fcy (vorticity) for a Cyclone case 

how the maximum parallel gradient, that is generated dynamically, depends on 
the parameters. The interest of the field- alignment technique would be some- 
what reduced if one were to discover, for example, that, for a given class of 
models, the maximum parallel gradient increases when the scale separation pa- 
rameter p* decreases. One must also stress that this behaviour may be model 
dependent. In particular, in the case of fluid models, it may depend on the par- 
allel closure. A Landau closure as the one employed here is somewhat weaker 
that a diffusive-type closure and it could allow a more frequent development of 
fine structures in the parallel direction. It would also be extremely interesting 
to study the behaviour of a kinetic model in this respect. 

6. Conclusions 

The main results of this work are here summarised. Field aligned techniques 
are necessary to optimise plasma turbulence codes. This work has shown that 
there are only two practical ways to achieve this. In the one employed by most 
codes, the position along the field lines is labelled by the poloidal angle like in 
the ballooning approximation of linear theory. In the one pursued here, points 
along field lines are labelled by the toroidal angle. 

This results in a more natural discretisation of the fields that has several 
advantages. In particular, it can deal with X-point geometries and both closed 
and open field lines. Furthermore, as shown in the Appendix, the description 
in the poloidal plane, where perpendicular derivatives are computed in most 
situations, need not employ magnetic coordinates. Information on field lines is 
needed only to compute parallel derivatives. This opens the way to generic and 
efficient coding for global tokamak simulations. 
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It has also been shown how conventional codes that employ a spectral, uni- 
form, representation on a magnetic surface can be converted to the new coordi- 
nate system with minimal programming. This has been done for the ETAI3D 
code, which now exists in two versions that solve exactly the same model equa- 
tions without further approximations. 

The comparison between the two versions has uncovered an unexpected prop- 
erty of the plasma turbulence model used for testing purposes. It has indeed 
been shown that parallel gradients can be intermittently higher than expected 
in the nonlinear phase. Parallel gradients are limited by the parallel heat flux, 
which, in this model, takes the form of a Landau closure. Landau closures are 
often used by gyro-fluid models. Then the key question, which cannot be an- 
swered here, is whether intermittently high parallel gradients are a real property 
of the weakly collisional plasma turbulence or an indication that the closure on 
the parallel heat flux is not adequate enough. Obviously, the whole construct 
of field- alignment rely on uniformly small parallel gradients for an efficient im- 
plementation. 
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Appendix A. Appendix 

This appendix sketches how the approach to field-aligned coordinates de- 
scribed in this work can be extended to avoid the use of magnetic surface coor- 
dinates to discretise the fields. The fields can then be discretised on a given grid 
related to the laboratory (or machine) reference frame. For a tokamak, these 
are the usual (R, Z) cylindrical coordinates such that Z is the direction of the 
torus axis, R the distance from the axis and $ the toroidal angle. The method 
is built on two concepts developed here: that micro-turbulent fields need to be 
known with high resolution only in two directions, and that knowledge of the 
magnetic field structure is needed only to compute parallel derivatives, which 
can be obtained by finite differences such that values at end points are obtained 
by interpolation. 

In the following, one considers a simple static, low-/3, cylindrical equilibrium, 
such that the suitably normalised magnetic field is given by 

B = b(x)+z (A.l) 

where one employs a Cartesian reference system (x, y, z) such that z is the 
direction of the magnetic axis, the main magnetic field along z is constant and 
normalised to unity, and b(x) is the poloidal magnetic field in the poloidal plane 
(x, y). The vector x indicates the position in this plane. The poloidal field can 
be written in terms of a flux function ^>(x) such that 

b = V x (ijiz) (A.2) 

Magnetic surfaces can be labelled by the value of ip. Both closed and open field 
lines can be treated. The parallel derivative operator is given by 

V|| = b V + d/dz (A.3) 

One has to look for a change of coordinates from the original (x, y, z) to a new 
set (£ a , s) such that s can be treated as a slowly- varying coordinate and only the 
two £ Q s [a = 1, 2) carry the information on the small scales. Taking advantage 
from what was learnt in Sec. [3j one considers a transformation of the form: 

C = l/ Q (x) + C"(x)(z-z fe ) (A.4) 
s = z — Zk (A. 5) 

where k labels a given sector in the z direction and V"(x) and C a (x) are yet 
unknown functions. In terms of the new variables the parallel derivative is given 

*>-» ! £w +i -* v '£* +o '* + k (A - 6) 
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In order to eliminate the fast-varying derivatives one has to satisfy the condi- 
tions: 



8V a 

° a = ^ (A-7) 

b a ^—=0 (A.8) 

dx a v ; 

In the following, the Poisson bracket notation is used such that the operation 
b a -§^k = [ip,A] for any function A. Eqs. (|A.7IA.8|) can then be written as 

[V,[^m=0 (A.9) 

Consider a function x( x i ) sucn that [ip 7 x\ — 1> whose solution can be found 
with the method of characteristics. This function identifies the position on 
ip = const surfaces and plays the role of a poloidal angle. Then a general 
solution of Eq. IA.9I is 

V a = r(^)+g a Wx(x 7 y), (A.10) 

where f a and g a are arbitrary functions. Thus, suitable solutions to Eqs. (|A.7I 
IA.8|) are found, such that the parallel gradient reduces to 

d 

v » = 

computed at fixed £ Q . 

In an actual code implementation based on finite differences, one defines any 
field at nodes in the Cartesian (x,y,z) grid. Derivatives in the poloidal plane 
are computed in this reference frame by holding z constant. 

In order to compute the parallel derivative by finite differences, one has 
to use function values at points (x + Ax, z% + Az) corresponding to a given 
increment As along s. This means finding end points along field lines for a given 
displacement Az along z. From Eqs. (|A.4HA.5|) one finds the finite difference 
equations for the unknown increments Ax 

[t{i>) + <?" Wx(x)] x+ Ax - [ 5 Q (V)]xAz = TO) + 5 " Wx(x)] x (A.ll) 

Solutions to these equations exist such as ^(x + Ax) = and x( x + ^ x ) = 
x(x) + Az. Thus end points for FD computations are obtained by moving along 
field lines for a given increment along z. 

Finally, in the special case of circular flux surfaces, such that ip = ip(x 2 +y 2 ), 
one finds that x is proportional to the usual poloidal angle, via the safety factor. 
Then the scheme sketched here reduces to the one described in Sec. [3] 
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